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We describe and extend the formalism of state-specific analytic density matrix renor¬ 
malization group (DMRG) energy gradients, first used by Liu et al (J. Chem. Theor. 
Comput. 9, 4462 (2013)). We introduce a DMRG wavefunction maximum overlap 
following technique to facilitate state-specific DMRG excited state optimization. Us¬ 
ing DMRG configuration interaction (DMRG-CI) gradients we relax the low-lying 
singlet states of a series of trans-polyenes up to C 20 H 22 . Using the relaxed excited 
state geometries as well as correlation functions, we elucidate the exciton, soliton, 
and bimagnon (“single-fission”) character of the excited states, and find evidence for 
a planar conical intersection. 


I. INTRODUCTION 

The density matrix renormalization group 1 ', introduced by WhittP, has made large ac¬ 
tive space multireference quantum chemistry calculations routine. In the chemistry context, 
there have been many improvements to White’s algorithm, including orbital optimization®*®, 
spin-adaptation- 2 ®, dynamic correlation treatments 1 ®®, and response theories 1 -^®, to name 
a few. The DMRG has been applied in many different electronic structure problems, ranging 
from benchmark exact solutions of the Schrodinger equation for small molecule d 16 * 1 ® ®, to 
active space studies of conjugated 7r-electron systems 2 ®®, the elucidation of the ground- and 
excited states of multi-center transition-metal clusterd®®, computation of high-order corre¬ 
lation contributions to the binding energy of molecular crystald®, relativistic calculations®, 
and the study of curve crossings in photochemistry®. 

Energy gradients are crucial to electronic structure as they define equilibrium structures, 
transition states, and reaction trajectories. Analytic energy gradients, introduced by Pu- 
la-ytHMU, are p re f erre( q numerical gradients, due to their low cost and numerical stability, 
and are now implemented for many single- and multi-reference quantum chemistry meth- 
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odJ™. 

Analytic DMRG energy gradients were first used by Liu et al. 36 in a study of the pho- 
tochromic ring opening of spiropyran. The theory, although simple, was not fully discussed. 
A contribution of the current work is to provide a more complete exposition of the theory 
behind DMRG gradients. A second contribution is to discuss their practical implemention 
in excited state geometry optimization. The simplest formulation of the gradients arises 
when the excited states are treated in a state-specific manner (that is, without orthogonal¬ 
ity constraints to lower states). However, such DMRG calculations can be susceptible to 
root-flipping, for example, near conical intersections. Furthermore, DMRG wavefunctions 
are specified both by a choice of active space as well as by discrete sets of quantum numbers 
associated with each orbital (used to enforce global symmetries, such as the total particle 
number). During an energy minimization, it is important that the wavefunction changes 
smoothly. Here, we present a state-specific DMRG analytic gradient algorithm that uses a 
maximum overlap technique both to stably converge excited states, and to ensure adiabatic 
changes of both the orbitals and the DMRG wavefunction during geometry changes. 

Trans-polyenes are well-known examples of molecules with interesting ground- and excited 
state structure, and form the central motifs for a large set of of biological compounds, such 
as retinals and the carotenoids. Many calculations using semi-empirical models (such as the 
Pariser-Parr-Pople (PPP) model) as well as with ab-initio methods such as multi-reference 
self-consistent-field (MC-SCF), have been carried out to identify the low-lying electronic and 
geometric features 48 ''. These studies, in general, give the following qualitative picture for 
long even polyenes: 1) Excitations, coupled to lattice relaxation, break the dimerization of 
the ground state and lead to local geometrical defect^ 78 . For example, adiabatic relaxation 
gives rise to a central polaronic feature in the bright 1 4 B U stat d 65 * 69 * 75 * 76 ( as well as two-soliton 
or four-soliton structures for the dark 2 4 A g stat^ 4 5665 . 2) For short polyenes ranging from 
ethylene to octatriene, studies of low-lying excited state relaxation pathways suggest that 
non-planar molecular conformations are important at energy crossings 79 95 . This opens up 
the question of the nature of energy crossings and their associated geometries in the excited 
states of longer polyenes. Here, using an ab-initio Hamiltonian and DMRG-CI analytic 
energy gradients, we revisit these questions in the excited states of relatively long trans- 
polyacetylenes, aiming for a more quantitative picture. 


In Sec. [TT] and HI we start by reviewing analytic energy gradient theory and DMRG 
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theory. In Sec. IV we discuss the formulation of DMRG analytic gradients. In Sec. |VJ we 
discuss the DMRG maximum overlap method for stable state-specific excited state calcula¬ 


tions without orthogonality constraints. In Sec. VI we apply our DMRG gradient algorithm 
to characterize the geometry and nature of the low-lying singlet states of the trans-polyenes. 


II. GENERAL ENERGY GRADIENT THEORY FOR VARIATIONAL METHODS 


For completeness, we briefly recall analytic energy gradient theory for variational wave- 
functions. The energy referred to is the Born-Oppenheimer potential, the sum of the elec¬ 
tronic energy and the nuclear-nuclear repulsion. The nuclear-nuclear repulsion gradient is 
trivial. The electronic energy (T|//|T) is the expectation value of the electronic Hamiltonian 

H, 


H 'y ^ tjj c ia C jo ~f y ^ Vijkl c ia c jcr/ c ka/ c lai 


( 1 ) 


ija 


ijklcrcj! 


where i, j, k,l are orthogonal (for example, molecular) orbital indices, and a, a' = {t, !}• H 
depends on the orbital functions through the integrals tij and Vijki- In a typical implemen¬ 
tation, these orthogonal orbitals are represented as a linear combination of atomic orbitals 
(AO’s) with orbital coefficients C: \i) = C ifl |/r), where we use Greek letters (/r, u) to 
denote AO orbitals. The geometry dependence of the integrals arises from the AO func¬ 
tions, which (as Gaussian type basis functions) explicitly depend on the nuclear positions, 
as well as through the LCAO coefficient matrix C, which also changes with geometry. The 
Hamiltonian may be regarded as a function of the nuclear coordinates {a*} and the coef¬ 
ficient matrix C: C). As the electronic wavefunction T) depends on variational 

parameters {c*}, the electronic energy is a function of {a*}, C, and {cj}: E({di}, C, {c,}). 
The energy gradient with respect to nuclear coordinate a takes the form 

dE (hi 


dE 

da 


dE 

da 


dE dC ifl 


Ifl 


dCi M da 




dci da 


( 2 ) 


If |T) is determined variationally, the energy is stationary to changes of {q}, thus the third 
term vanishes. The gradient then only depends on the change in nuclear coordinates and 
orbital coefficients, 


dE 

da 


dE 
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dE dCi 
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( 3 ) 
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It is convenient to rewrite the energy gradient in terms of density matrices. The energy 
is expressed as 

E ^ ' t jj'Jjj T ^ ^ Ujjkl^ijkl i (4) 

ij ijkl 

where 7 ^ = (T c\ a c ja |T) and Y ijkl = \ (^1 c L c ]ai c kai c ia I'h) are the one- and two- 

particle density matrices. Since the second-quantized operators have no dependence on 
geometry, and the wavefunction depends only on {c*}, it follows from Eq. ([3]) that the 
energy gradient is expressed in terms of the one- and two-electron derivative integrals and 
density matrices, 

dE \ dtj/j x (Ivj jki ,, 

(5) 


dE 

da 


E dtij \ dvij k i 

~^ lij + 2^ ~^r Tijkl - 


ij ijkl 

The one- and two-electron derivative integrals involve the orbital derivative (response), 
dC/da. Writing 


da 


l/l 


da 




( 6 ) 


gives 


dE 

da 


E dhij dv^ki 


E x ‘r S,i 


i ijkl 

E - vo. 


da 
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(7) 


where S t j = ( i\j ) is the overlap matrix of the orthogonal orbitals i and j, 

Sij = ''j ^ CifjS^pCjy, 

/111 


( 8 ) 


and S /w = (/i \v) is the overlap matrix in the underlying AO basis. X is the Lagrangian 
matrix in the Generalized Brillouin Theorem (GBT)®, given as 


Aj j ^ ^ hj jn'lmj T 2 ^ ^ i-Urnkl Yjrnkl i 


(9) 


mkl 


characterizing the energy cost of electronic excitation from itli orbital to jth orbital. The 
gradient formula can be rewritten entirely in terms of AO quantities 4 -^, 


dE \ dhn V \ 

£>-wr + £ r 
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where h /JM , 7 ^, v pvpa , Tp Vpu are the one- and two- particle integrals and reduced density 
matrices, respectively, in the AO basis. 

In general, the orbital derivative dCi^/da requires the solution of equations which couple 
the wavefunction coefficients <7 to the orbital coefficients C t/l . However, there are two com¬ 
mon situations where the orbital response is simplified. The first is when the orbitals are 
defined independently of the correlated wavefunction, for example, for Hartree-Fock (HF) 
or Kohn-Sham (KS) orbitals. Using the Hartree-Fock canonical orbitals as an example, the 
orbital response U a is defined by the Hartree-Fock convergence condition, 


dF i:j 

da 


0 


( 11 ) 


which leads to the definition of the U a matrix 


vir d.o. 


u a ■ = 

13 


YYxvi ; +ay, 


i c v 


( 12 ) 


where 


■A-ij^kl ^Vijkl Vikjl Vnjk 

s.“- = F?j- Stfy - £ S m(- v w)’ 

jk 


(13) 


with F°- = dFij/da, = dSij/da, and the various e are HF orbital energies. Eq. (12) is 


the coupled-perturbed Hartree-Fock (CPHF) equation®', and uniquely defines the U a matrix 
elements for canonical orbitals. In a similar way, other types of orbital response, for example 
for the Kohn-Sham orbitals, or localized Hartree-Fock orbitals, can be computed from the 
corresponding coupled-perturbed single-particle equation ^ 41 98 . 

The second simplifying case is when the correlated wavefunction energy is itself stationary 
with respect to orbital variations. In this case X t] — Xji = 0, and the orbital response is 
not required, even though it is formally coupled to the correlated variational wavefunction 
coefficients. The energy gradient reduces to the simpler form, 
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III. GENERAL DMRG THEORY 

The DMRG is a variational wavefunction method". For a set of L orthogonal orbitals 
(where the states for the ith orbital are |cr,) = {| 0 ), |t) , |i) , we choose a partitioning of 

the orbitals into a left block, single site, and right block, consisting of orbitals {1...Z — 1 }, {/} 
and {/ + respectively. The corresponding canonical “one-site” DMRG wavefunction 

takes the matrix product form 

|T) = L ai L rT 2 ...L a ‘- 1 

0 - 10-2 ...a L ( 15 ) 

x C ai R ff|+1 R cr,+2 .. .R 0 ^ l ) . 

The (rotation) matrices L CT,: and R CTi are of dimension M x M, except for the first and 
last which are of dimension 1 x M and M x 1 respectively. They satisfy the left- and 
right-canonical conditions 

L aiT L ai = i 

Y R CTi R CT * T = 1 (16) 

while the C a ‘ (wavefunction) matrix satisfies the normalization condition 

tr^C CTiT C ffi = 1. (17) 

O'Z 

Together, {C CTi } and {R CTi } contain the variational parameters. As in other variational 

methods, the coefficients of the matrices are determined by minimizing the energy. In 
principle, a direct gradient minimization of the energy with respect to all the matrices, 
subject to the canonical conditions Eq. may be performed. In practice, the DMRG 

sweep algorithm is normally used. Here, at a given step l of the sweep, corresponding to the 
block partitioning {1...1 — 1}, {/} and {/ + the energy is minimized only with respect 

to the C CTi wavefunction matrix, with the {Id 7 *}, {R*} rotation matrices held fixed. The 
minimizing C ai is obtained from an effective ground-state eigenvalue problem 

He = Ec (18) 

where c denotes C ai flattened into a single vector, and H denotes H expressed in the basis 
of renormalized basis states defined by the {L^}, {R^} matrice^. In the next step of the 
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sweep, the single site is moved from l to l +1 (or l to l — 1 in a backwards sweep). To satisfy 
the new canonical form with the single site at l + 1, where the C ai matrix is replaced by 
an L a ‘ matrix, and the l + 1 site is associated with a new C ai+1 matrix, we use the gauge 
relations, 


C ai = L CTi A 

C ai+1 = ART+ 1 . (19) 

By sweeping through all the partitions l — 1... L, and minimizing with respect to the C ai 
matrix at each partition, the DMRG sweep algorithm ensures that all the variational degrees 
of freedom in the DMRG wavefunction are optimized. 

An important aspect of DMRG calculations is the enforcement of symmetries, including 
global symmetries such as the total particle number and spin. In the DMRG wavefunction, 
Abelian global symmetries (such as total particle number) are enforced by local quantum 
numbers. For example, to enforce a total particle number of N in the wavefunction, each 
value of the 3 indices a, i,j in the matrix elements L?-, R" , C"- can be associated with 
an additional integer N i} Nj, N a . (These values can be interpreted in terms of the particle 
numbers of the renormalized states (for N t and Nj) and for the states of the single site (for 
N a )). Then, a total particle number of N is enforced with the rules: 

L : Ni + N a = Nj 
R :Nj + N a = N 

C:Ni + N a + Nj = N. (20) 

Applying these conditions to L?■, R?■, C^- means that the matrices have a block-sparse 
structure, which is important to maintain during geometry optimization. 


IV. STATE-SPECIFIC DMRG ANALYTIC ENERGY GRADIENTS 

At convergence of the above (one-site) DMRG sweep algorithm, the contribution of the 
wavefunction coefficients to the gradient ( dci/da in Eq. (|2]) ) vanishes, as expected for a 
variational wavefunction method. Thus the analytic energy gradient theory for variational 
wavefunctions described in Sec. [TTjcan be applied. 


We will consider energy gradients for two kinds of DMRG calculations. The first are 
DMRG configuration interaction (DMRG-CI) analytic gradients, using HF canonical or¬ 
bitals. In this case, the orbital response is given by the CPHF equations, presented in 
Sec. [TTJ The DMRG calculations are carried out within an active space, chosen as a subset 
of the canonical orbitals. Because the DMRG wavefunction is not invariant to rotations of 
the active space orbitals for small M, the contribution of the active orbital response must 
be computed specifying a particular orbital choice (rather than just their manifold), such as 
the canonical HF orbitals. 

The algorithm to compute the DMRG-CI analytic gradient with HF canonical orbitals is 
as follows: 


1. Solve the HF equations for the canonical orbital coefficient matrix C. 

2 . Select an active space, and solve for the DMRG wavefunction in this space. Compute 
the one- and two-particle reduced density matrices 7 \j and at the convergence of 
the single-site sweep algorithm. 


3. Compute the AO derivative integrals dh^/da, dv pupa /da and dS pu /da, and the X 
matrix in Eq. 0 - 


4. Use the derivative integrals to construct the CPHF equation in Eq. (12) (or the equiv¬ 


alent Z-vector equation 42 ) and solve for U a for all nuclear coordinates. 


5. Compute the energy gradient by contractions of all the above integrals and matrices 
according to Eq. 0 or gof. 


The second kind of DMRG calculation we consider is a DMRG complete active space 
self-consistent field (DMRG-CASSCF) calculation. For DMRG-CASSCF wavefunctions, 
the DMRG energy is stationary to any orbital rotation, thus 


A 'ij — Xji — 0 


( 21 ) 


and by Eq. 0 and (10) this means that the orbital response is not required even though it 
is coupled to the response of the DMRG wavefunction. However, because the DMRG wave- 
function is not invariant to active space rotations for small M, it is necessary to optimize 
the active-active rotations also, unlike in a traditional CASSCF calculation. Alternatively, 
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if active-active rotations are omitted, the DMRG-CASSCF gradient can be viewed as an ap¬ 
proximate gradient with a controllable error from active-active contributions (which vanishes 
as M is extrapolated to oo.) 

The algorithm for the DMRG-CASSCF gradient is: 

1. Solve for the DMRG-CASSCF orbitals with the one-site DMRG wavefunction. In each 
macroiteration: 

(a) Solve for the one-site state-specific DMRG wavefunction, and compute the one- 
and two-body reduced density matrices 7 ^ and r^. 

(b) Using 7 ij and T^i, compute the orbital gradient and Hessian, both of which 
include elements for active-active rotations. 

(c) Update the orbitals with the orbital rotation matrix. 

2 . Compute the AO density matrices 7 ^ and V pvpa at the convergence of DMRG- 
CASSCF. 


3. 

4. 


Compute the AO derivative integrals dh pv /dA, d(pv\pa) / dA and dS pv /dA. 


Contract all the above integrals and matrices using Eq. (14) to obtain the energy 
gradients. 


V. ADIABATIC ORBITAL AND WAVEFUNCTION PROPAGATION AND EXCITED 
STATE TRACKING 


Geometry optimization requires adiabatically propagating along a potential energy sur¬ 
face. For a DMRG calculation, this means that in each geometry step, the orbitals defining 
the active space should change continuously, and the quantum numbers and associated 
block-sparsity pattern of the matrices should not change. The former can be achieved using 
maximum overlap techniques, while the latter can be done by fixing the quantum numbers 
at the initial geometry. For state-specific excited state calculations, the maximum over¬ 
lap technique is further important to prevent root-flipping. Root flipping in state-specific 
DMRG calculations arises because the matrices optimized in the wavefunction for one state 
(15) are not optimal for another staW 100 . (Note that the gradient formalism presented above 


is only valid for state-specific, rather than state-averaged, DMRG calculations). 
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A. Orbital maximum overlap 


The maximum overlap technique for the orbitals involves computing the overlap matrix 
between MO’s of the (m — l)th and mth step 


rim—l.m / / 

S<4 = W 


m— 1 1 n i,m 


13 


m 


= 2 .^ 


l CTA4> r'lO 


( 22 ) 


where (</>™ are the AO overlap matrix elements of (m — l)th and mth steps. For the 

active space, we choose the orbitals at step m with maximum overlap with the active space 


orbitals at step m — 1. Eq. (22) also allows us to align the MO phases for adjacent geometry 
optimization steps. 


B. Excited state tracking in DMRG 

We further use maximum overlap of the DMRG wavefunctions to target and track the cor¬ 
rect state-specific excited state solution. Within the standard ground-state sweep algorithm 
at a given geometry, the desired excited state can usually be found in the eigenspectrum at 
the middle of the sweep (when the renormalized Hilbert space is largest) but can be lost at 
the edges of the sweep when the renormalized Hilbert space is small (if it is generated for 
the incorrect eigenvector). To keep following the excited state across the sweep by generat¬ 
ing the appropriate renormalized Hilbert space, we ensure that at each block iteration we 
always pick the Davidson solution with maximum overlap with the excited state solution at 
the previous block iteration. Between geometries, we ensure that we are tracking the correct 
excited state by computing the overlap between the DMRG wavefunctions at the different 
geometries. In principle, this requires multiplying the overlaps between the L CT , R CT matrices, 
and c vectors. However, we find it is sufficient (and of course cheaper) to only compute the 
overlap between the c vectors for the two geometries, at the middle of the sweeps. 

The state-specific DMRG wavefunction maximum overlap scheme is: 

1 . At the initial geometry, use a state-averaged DMRG algorithm to obtain initial guesses 
for n stated 00 . (The more robust two-site DMRG algorithm may be used here 1 ®, and 
a highly accurate initial guess for a small M can be obtained by running back sweeps 
from large M 102 ). Store the wavefunction vectors {c*} (for i — 1,2,..., n) at the middle 



11 


of the sweep. Note that in the state-averaged procedure all n states share the same 
left and right rotation matrices {L 17 } and {R CT }. 

2 . At a new geometry optimization step (=initial geometry in the first step), restart 
the DMRG sweep with the same M from the previous solution for the targeted ex¬ 
cited state, and use state-specific DMRG with the one-site sweep algorithm to get 
the new solution for the targeted excited state. (Note that any noise in the DMRG 
algorithm should be turned off). At each block iteration, apply the following steps in 
the Davidson solver: 


(a) Perform DMRG wavefunction prediction by Eq. (19) from the previous block 
iteration, to obtain guess vectors {c l ss } for the current block iteration. 


(b) Perform the Block-Davidson algorithm to obtain solutions {c* oi }. 

(c) Compute overlaps between vectors {c* oZ } and {c l guess }, and align the phases when 
needed. 


(d) Choose the new solution c x sol in {c l sol } for the targeted excited state, from the 
largest overlap between c x sol and c™ uess . 

(e) Store the vector c x sol as the new solution. 


3. Repeat Step. [2] in further geometry optimization steps. 


VI. EXCITED STATE GEOMETRIES OF TRANS-POLYENES 

Excited state geometry optimization in linear polyenes serves as a starting point to un¬ 
derstand the photophysical and photochemical behaviour of analogous systems, such as 
the carotenoids, in biological processes. We take as our systems, the trans-p olyacetylcnes 
C 2 nH 2 n+ 2 , with n — 5 — 10. We modeled the excited states and geometry relaxation as fol¬ 
lows: 1) We obtained ground state S 0 (1 1 A 9 ) geometries with DFT/B3LYFG° 3j . 2) We then 
used the DFT ground state geometries as initial guesses to perform ground state geometry 
optimization with DMRG-CI analytic energy gradients. 3) We recomputed excited states at 
the DMRG optimized ground state geometries. 4) We then further relaxed the excited state 
geometries with the DMRG-CI gradients. All calculations were performed with the cc-pVDZ 
basis sell® 4106 . The active spaces were chosen as (ne, no), where n is the total number of n 
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electrons. We identified the i r active spaces consisting of carbon 2 p z orbitals from the Lowdin 
MO population analysis at the initial geometry, and tracked the active spaces through the 
geometry relaxation with the orbital maximum overlap method in Sec. V A We also car¬ 
ried out additional calculations with a second “energy-ordered” active space, consisting of 
the lowest it and a orbitals to make up an (ne, no) active space. We clearly distinguish 
when we are referring to the second active space in the discussion below. The initial ground 
state DFT/B3LYP geometry optimizations were carried out with the Molpro quantum 
chemistry package 107 . State-specific DMRG wavefunctions were obtained with the Block 
DMRG program® 13 S3 using the state-specific and adiabatic wavefunction tracking by 


wavefunction maximum overlap in Sec. |VB| DMRG-C1 gradients were implemented in the 
ORCA quantum chemistry package. All calculations worked in the canonical HF orbital 
basis (no localization). To improve the geometry optimization we employed approximate 
nuclear Hessians, updated by the BFGS method 109112 . 

To simplify the analysis, in this work we only considered geometry optimization in the 
plane. Non-planar geometries are of course relevant to polyene excited states but even at 
planar geometries, important features of the electronic excited state geometries (e.g. the 
solitonic structure) appear and remain to be understood at an ab-initio level. The planar 
optimization was not enforced explicitly other than through a planar initial guess, and other¬ 
wise the coordinates were allowed to relax in all degrees of freedom. Consequently, electronic 
wavefunctions were computed within C\ spatial point group symmetry. We used three dif¬ 
ferent numbers of renormalized states M=100, 500, 1000 to obtain DMRG wavefunctions for 
all states, to examine the influence of wavefunction accuracy on the geometries. Converg¬ 
ing DMRG wavefunctions to a high accuracy ensures the accuracy of the particle density 
matrices, which then ensures that the correct geometric minima can be reached. However, 
when the magnitude of gradients was much larger than the unconverged DMRG error, (for 
example, when the geometry was far from the equilibrium) loose DMRG convergence and 
fewer sweeps were used to decrease the computational time. 

To further characterize the low-lying excited states, we analyzed the exciton and bi- 
magnon character of state transitions using the transition particle density matrices. The 
first-order transition density matrix element in the MO basis between the ground (GS) and 
excited states (ES) is 




ES 


c\ c j 




GS 


(23) 
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where i and j denote spatial MO indices. We used the first-order transition density matrix 
to locate the first optically dark and bright states by the following well established state 
signatures: 1) A single large element where ?'=LUMO, j=HOMO, indicating the first op¬ 
tically bright state. 2) Two dominant elements where i = LUMO + 1, j — HOMO and i 
= LUMO, j = HOMO - 1 indicating the first optically dark state. Real space particle-hole 
excitation patterns were further analyzed by the real space first-order transition density 
matrix, which was obtained by transforming the vir-occ block of the MO first-order transi¬ 
tion density matrix to the orthogonal 2 p z basis. Real space particle-hole excitation patterns 
were characterized by excitations of an electron from an orbital at R — r /2 to an orbital 
to R + r/2, where R was set at the centre of a polyene chain, and r is the particle-hole 
separation length. We illustrate the excitons graphically by plotting (cj,c n _ p ), where p is 
the index of the carbon 2 p z orbital, and n is the total number of 2 p z orbitals in the chain. 

Similarly, the real space bimagnon character is characterized by the real-space double-spin 
flip transition density 

ES Cp’O-Cp'—o-Cn—p —uCn—p^cr ^ GS ^ (24) 

where a = {t, i}- The real space second-order transition density matrix was transformed 
from the vir-vir-occ-occ block of the MO basis second-order transition density matrix. 

Analogously to previous studies, we further examined bond orders and geometrical defects 
(solitons) through the bond length alternation (BLA) function 5 n 

S n = (- l) n+ \x n+1 -x n ) (25) 

where n = 1, ...,N bond , and x denotes bond lengths. For even-site trans-polyacetylenes, the 
two edge bonds at the ground state are always double bonds, thus 5 n will always be positive. 
Consequently, negative values of S n indicate a reversed bond order, and a vanishing ( S n = 0) 
value comes from two equal bond lengths, i.e., an undimerized region. 

A. State signatures and geometries 
1. Ground state So 

The ground state of polyenes is denoted by the symmetry label 1 1 A g (here we are using 
symmetry labels characteristic of idealized symmetry) and the relaxed ground state 
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FIG. 1. Bond length alternation function 5 n for relaxed ground state geometries of C 20 H 22 , from 
left to right. 


geometries are planar and dimerized. For the ground state, DMRG wavefunctions with M 
= 100 are sufficient to achieve qualitative accuracy in bond lengths of our studied polyenes. 
M = 500 is sufficient for quantitative accuracy. For example, M = 100 produced errors of 
no more than 0.006 A for C 20 H 22 , while M = 500 converged the bond lengths to an error 
of 0.0003 A, as compared to bond lengths using M = 1000 (near exact). This finding is 
consistent with the ground state wavefunction of even-carbon trans-polyenes being mostly 
a single-determinant, and thus accurately described by DMRG in the canonical molecular 
orbital basis with small M. 

The BLA function S n of the relaxed ground state geometry of C 20 H 22 from DMRG and 
DFT is shown in Fig. [l] The BLA functions from both DMRG and DFT give the same 
pattern, showing a weaker dimerization in the middle region compared to the edges of 
the carbon chain. Compared to the dimerization in DFT, the dimerization in the DMRG 
calculations is suppressed, indicating a smaller dimerization gap. Compared to DFT, a Tr¬ 
active space Hamiltonian (as used in the DMRG calculations) is associated with a larger 
effective Coulumb interaction U due to the lack of dynamic correlation. A suppression of 
the dimerization can then be expected, as the dimerization magnitude behaves as f/~ 3 / 2 in 
the strongly interacting limit 113 . 
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FIG. 2. Vertical and relaxed excitation energies from DMRG optimized geometries: vertical So-Si 
(opened squares), vertical S 0 -S 3 (opened circles), relaxed S 0 -S 1 (solid squares), relaxed S 0 -S 3 (solid 
circles). 

2. Excited states 


The first optically bright state in the single -tt complete active space is the third excited 
state S 3 , and denoted by the symmetry label 2 1 B U . The corresponding MO based first-order 


transition density matrix between S 3 and So (defined by Eq. (23)) possesses an element ~ 1.0, 
where i = LUMO and j = HOMO, along with other elements < 0.1. This signals a (HOMO 
—y LUMO) single particle-hole transition, characteristic of the first optical transition. 

The 1 1 B„ state corresponds to the second excited state S 2 in the single-7r active space. 
Notable first-order excitations in the So/S 2 transition, for instance in CioH 12 , are (HOMO —y 
LUMO + 2) and (HOMO - 2 —* LUMO) excitations, both with elements ~ 0.5 at the ground 
state equilibrium geometry. A large (HOMO —> LUMO) excitation is missing for the So/S 2 
transition for all the polyenes, ruling it out as the usual bright state. Note that the order 
of excited states depends on the choice of active space, i.e., the effective Hamiltonian. If 
one changes from the single- tt active space to an energy-ordered active space which includes 
both cr and 7 r orbitals within the (ne, no) active space window, one finds that the 1 1 B U state 
is an S 2 state corresponding to the physical optically bright HOMO —> LUMO transition. 
This demonstrates the well-known strong effect of dynamical correlation on the low-lying 
excited state order in linear polyenes. 

The first optically dark state is the Si state, denoted by 2 1 A g . The S 0 /Si transition 
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exhibits dominant (HOMO —> LUMO + 1) and (HOMO - 1 —> LUMO) single excitations, 
along with a dominant (HOMO, HOMO—» LUMO, LUMO) double excitation. The position 
of this low-lying excited state remains as the Si in an energy-ordered active space. 

For the bright state, optimized bond lengths were not strongly dependent on M. For a 
small system such as Ci 0 H 12 , M = 100 produced a largest error of 0.0003 A in the bond 
lengths, as compared to the M = 1000 result. For a larger system such as C 2 oH 2 2, the 
bond lengths at M = 100 differed the ones at M = 1000 by no more than 0.005 A, and 
the largest error at M = 500 was only 0.0006 A. For the dark state, however, the precision 
of the optimized geometry was more sensitive to the choice of M for the longer polyenes. 
This may not be surprising, as the first optically bright state is mainly a single-reference 
state, while the lower dark state has more challenging multi-reference character 114 . For all 
the polyenes considered, if we use small M, the largest error in the bond lengths of the dark 
state occurs for bonds around the geometrical defects (solitons). In C 2 oH 22 , the largest error 
at M = 100 is about 0.025 A, coming from the bonds C3-C4 and C16-C17 which are around 


the solitons (see in Sec. VIC). On the other hand, central bonds in the dark state are much 
less dependent on M, e.g. M = 100 yields errors < 0.012 A for bonds from C 6 -C7 to C13-C14 
in C 2 oH 22 . In a localized real space view, this behaviour reflects the strong localization of 
multi-reference electronic structure around the geometrical defects. 


B. Excitation energy 

We show vertical and relaxed excitation energies as a function of 1/n for the first optically 
dark (2 1 A £ ,) and first optically bright (2 1 B U ) states for all considered C 2n H 2ri+2 in Fig. [2} 
Compared to the experimental excitation energies for CioHi 2 to Ci4Hi 6 in hydrocarbon 
solutiond 11 ^, our relaxed excitation energies are 0.3 eV higher for the relaxed dark state, and 
1.7 eV higher for the bright state. This is in part due to the lack of dynamic correlation in 
our calculations as well as basis and solvent effects. 

The dark state is always observed as below the bright state. We observe relaxation 
energies for all the polyenes of about 0.35 eV for the bright state and about 1.20 eV for the 
dark state. The substantial relaxation energy for the dark state is consistent with the much 
larger geometry relaxation as compared to the bright statd®2. 

Our calculations find the 1 1 B U state to lie relatively close to the 2 1 B U state at the ground 
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Geometry optimization step 


FIG. 3. Energies of the S 2 and S 3 states in the S 2 geometry optimization of C 20 H 22 , computed with 
the energy-ordered active space, as a function of the geometry relaxation step. At the ground state 
equilibrium (step 0), the S 2 and S 3 states are FEU (first bright state) and 2 1 B M states respectively. 
At step 4 the molecule gives a S 2 and S 3 gap of 0.019 eV, strongly indicative of a conical intersection. 
After this step the 1 1 B U and 2 1 B. U states are swapped in terms of the state energy order. (Note 
that the S 3 state energy oscillates as only the S 2 state energy is being minimized). The molecular 
geometry remains planar along the relaxation. 

state equilibrium geometry, with a l 1 B tt -2 1 B u energy gap consistently about 0.27 eV for all 
the polyenes. Given the small magnitude of this energy gap, it seems likely that there can be 
an energy crossing between 2 1 B U and 1 1 B U states. If we use the energy ordered active space 
we do find an energy crossing between these states for C 20 H 22 at a planar geometry, near the 
Franck-Condon region (Fig. [ 3 ]). Of course we also expect non-planar conical intersections, 
as previously found in butadiene 93 94 and octatriene 95 . 

C. Solitons 

The BLA S n functions for the first optically dark (2 4 A g ) and first optically bright (2 4 B U ) 
states are shown in Figs. [4] and [5j These curves are almost parallel across all the polyenes 
for the dark and bright state respectively, indicating generally similar behaviour across the 
systems. 

For the 2 3 A g state, the BLA in short polyenes C 10 H 12 and C 12 H 14 is completely reversed 
from the ground-state, as shown by the all negative 5 n values along the chain. The reversal 
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FIG. 4. Bond length alternation function 5 n for relaxed first dark state geometries, from edge (left) 
to center (right). 



FIG. 5. Bond length alternation function 5 n for relaxed first bright state geometries, from edge 
(left) to center (right). 

of BLA in 2 1 A g in short polyenes has previously been understood in terms of the dominant 
valence bond configurations 88 with reversed BLA. For long polyenes, undimerization emerges 
near the edges as shown by changes in the sign of the S n functions, and the BLA is opposite 
on the two sides of the undimerized regions. This result is in agreement with earlier semi- 
empirical studies on long polyenes 66 -®^, and our result shows the two-soliton structure in the 
relaxed 2 x A g state. 

For the 2 1 B U relaxed geometry, 8 n systematically shows a polaronic defect in the chain 
centre. This is also consistent with previous semi-empirical studies 66 67 . For short polyenes, 
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the vanishing dimerization in the central region can be understood in terms of ionic VB 
configurations along the chain! 58 ! In terms of excitons, the polaronic geometry is also viewed 
as evidence of a bound particle-hole excitation localized near the chain centre' 8 . 


D. Excitons 

Within the one-electron manifold, we can visualize the excitons with the real space 
particle-hole excitation density (c£c n _ p ). As we relax the geometry, we can observe the 
shape of the exciton change. Geometry relaxation is important to overcome the exciton self- 
trapping'^, e.g., in a polyene chain in its dimerized ground-state geometry. The real space 
particle-hole excitation densities of C 20 H 22 are shown in Fig. [6] and Fig [TJ for the bright and 
dark state respectively. 

At the ground state equilibrium geometry, i.e., a dimerized geometry, the particle-hole 
excitations of the bright state are strongly bound, as seen in Fig. [6|a). This is similar to as 
seen in the single-peak real space exciton structure from DFT-GWA-BSE calculation^ 116 ', as 
well as the n — 1 Mott-Wannier exciton pattern in the weak-coupling limit 66 . For the dark 
state, particle-hole pairs are slightly separated at the dimerized geometry, as illustrated 
by the double-peak real space exciton structure in Fig. |6^a). This has been identified in 
previous studie d 6 * 116 !, as an n = 2 Mott-Wannier exciton. However, the amplitudes of the 
densities are ten times smaller as compared to that of the bright state, essentially suggesting 
neglegiblc exciton character for the dark state reached by a vertical transition. 

After geometry relaxation, the particle-hole separation in the bright state increases, al¬ 
though the particle-hole pair remains bound at the bright state equilibrium geometry, as 
shown in Fig. |b]^b). For the dark state, however, geometry relaxation seems to unbind the 
particle-hole pair, as shown by largely separated peaks in Fig. [7^b). Along with the enhanced 
transition density amplitude, this suggests the emergence of a long-distance charge-transfer 
character associated with the dark state equilibrium geometry. 


E. Bimagnons and singlet fission in 2 1 A g 


The relaxed dark state 2 1 A £ , geometry possesses a separated two-soliton structure as dis¬ 


cussed in Sec. VIC The locally undimerized regions in the relaxed 2 X A g state can be thought 
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FIG. 6 . Real space particle-hole excitation density of C 20 H 22 between ground state and first bright 
state, computed at relaxed geometries of (a) ground state (b) first bright state. 


to arise from a form of “internal singlet fission 1 ' 7 ®, i.e., forming local triplets (magnons) while 
the total spin remains a singlet. The local triplets can be identifed from the local peaks of 
the real space spin-spin correlation function of the 2 l A g wavefunction as in Ref. 66 . 


Here, we can also characterize the bimagnon character by the real space double-spin flip 


transition density between the S 0 and 2 1 A g states (see Eq. (24)). We show the real space 
double-spin flip transition density of C 2 oH 2 2 as a function of the site index in Fig. [ 8 j At the 
ground state equilibrium, the bimagnons are confined near the chain centre, as indicated 
by the local central double peaks. However, the bimagnons are highly mobile, and with 
geometry relaxation, the singlet fission character becomes much more delocalized. The 
transition density distribution possesses two peaks at carbon 3 and 17 at the dark state 
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FIG. 7. Real space particle-hole excitation density of C 20 H 22 between ground state and first dark 
state, computed at relaxed geometries of (a) ground state (b) first dark state. 

equilibrium geometry, which is consistent with the positions of the solitons shown in Fig. |4} 


VII. CONCLUSIONS 

We presented the detailed formalism for state-specific DMRG analytic energy gradients, 
including a maximum overlap algorithm that facilitates state-specific excited state geom¬ 
etry optimizations. We employed these techniques to study the ground and excited state 
electronic and geometric structure of the polyenes at the level of DMRG-CI. Our quanti¬ 
tative results are consistent with earlier qualitative semi-empirical studies of the exciton, 
bimagnon, and soliton character of the excited states. In addition to complex bond-length 
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FIG. 8. Real space double-spin flip density between the ground state and the first dark state. 

alternation patterns, we find evidence for a planar conical intersection. 

DMRG analytic energy gradients provide a path towards the dynamical modeling of 
excited state and highly correlated quantum chemistry. The interaction of dynamic and 
non-adiabatic effects with strong electron correlation remains an open issue, which can now 
be explored with the further development of the techniques described here. 
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